home
***
CD-ROM
|
disk
|
FTP
|
other
***
search
/
Ham Radio 2000
/
Ham Radio 2000.iso
/
ham2000
/
satellit
/
usat92b2
/
usatqb.bas
< prev
next >
Wrap
BASIC Source File
|
1993-02-03
|
28KB
|
890 lines
' AMSAT ORBITAL PREDICTION PROGRAM de W3IWI - May 1980
' COPYRIGHT 1980 by Dr. Thomas A. Clark, W3IWI
' 6388 Guilford Road
' Clarksville, MD. 21029
' REVISED & MODIFIED FOR IBM-PC by R. D. Welch, W0SL - May 1983
' 908 Dutch Mill Dr.
' Manchester, Mo. 63011
' ENHANCED AND DEBUGGED BY Ing. H.F.STRASSER, OE1HSI - JAN. 1985
' A 1238
' VIENNA, AUSTRIA
' Modified to use the NASA/NORAD two-line element set format (as distributed
' by Dr. Kelso of the Air Force Institute of Technology) and revised for
' QBASIC by Antonio Querubin, AH6BW
' Internet: tony@mpg.phys.hawaii.edu
' AMPRnet: ah6bw@uhm.ampr.org
' PBBS: ah6bw@ah6bw.hi.usa.oc
' BITnet: querubin@uhunix
' Last update: January 15, 1993
' Permission granted for non-commercial use providing
' credit is given to the author(s), AMSAT and ORBIT Magazine.
DEFDBL A-Z
DECLARE FUNCTION gets$ ()
DECLARE SUB getkey (k$)
DECLARE SUB timezone (d0$, t0$, d1$, t1$, tzcor!)
DECLARE FUNCTION juliantodate$ (day%, year%)
DECLARE FUNCTION datetojulian% (d$)
DECLARE FUNCTION trim$ (i$)
DECLARE FUNCTION readtle% (tlefile$, name$(), epochyear%(), epochjulianday#(), inclination#(), epochRAAN#(), eccentricity#(), epochArgOfPerigee#(), epochMA#(), MeanMotion#(), epochrevolution&())
DECLARE SUB readgrnd ()
DECLARE SUB plotloc (longitude!, latitude!, xcoord%, ycoord%)
DECLARE FUNCTION arctan# (y#, x#)
DEF fnyear% = VAL(RIGHT$(d$, 2))
DEF fnhour% = VAL(LEFT$(t$, 2))
DEF fnmonth% = VAL(LEFT$(d$, 2))
DEF fnday% = VAL(MID$(d$, 4, 2))
DEF fnminute% = VAL(MID$(t$, 4, 2))
DEF fnsecond% = VAL(RIGHT$(t$, 2))
DEF fnt$ (d%) = MID$(STR$(d% + 100), 3)
DEF fnleap% (year%) = ((year% MOD 4) = 0)
DEF fnradians (degrees) = degrees / 180# * pi#
DEF fndegrees (radians) = radians / pi# * 180#
DEF fnarcsin (sine) = ATN(sine / SQR(1 - sine * sine))
' Greenwich Mean Sidereal Time
DEF fngmst# (YR%) = (99.6367# - .2387# * (YR% - 1989) + .9856# * INT((YR% - 1989) / 4#)) / 360!
CLEAR
maxsats% = 17
DIM name$(maxsats%)
DIM epochyear%(maxsats%), epochrevolution&(maxsats%), satl(maxsats%, 2)
DIM epochjulianday#(maxsats%), inclination#(maxsats%), epochRAAN#(maxsats%)
DIM eccentricity#(maxsats%), epochArgOfPerigee#(maxsats%)
DIM epochMA#(maxsats%), MeanMotion#(maxsats%)
DIM BeaconFrequency&(maxsats%), SemiMajor#(maxsats%)
DIM sat%(7), pkt%(7), kep%(7) ' Image arrays (9 x 5 pixels)
DIM cc#(3, 2) ' Translation matrix from the satellite orbital plane coordinate
' system to the celestial coordinate system (origin at geocenter)
' ****** NUMERIC CONSTANTS ******
pi# = 3.141592653589793# ' Value of PI
r0# = 6378.14# ' Earth's radius
f# = 1# / 298.16# ' 1/Earth flattening coef.
g0# = 75369793000000# ' GM of Earth in (orbits/day)^2/km^3
g1# = 1.0027379093# ' Solar/Sidereal time ratio
c# = 299792.458# ' Speed of light (km/sec).
' String constants
tlefile$ = "USAT.TLE" ' Default NASA/NORAD element file
' YOU WILL NEED THE FILES USATCGA.MAP, USAT.TLE AND USATGRND.DAT
' TO RUN THIS PROGRAM
' **** SATAUS Menuprogramm de OE1HSI VERSION 1.5 26.JAN. 1985
' ****** ESTABLISH SATELLITE ELEMENT MATRIX ******
PRINT "Reading satellite elements from "; tlefile$
numsats% = readtle%(tlefile$, name$(), epochyear%(), epochjulianday#(), inclination#(), epochRAAN#(), eccentricity#(), epochArgOfPerigee#(), epochMA#(), MeanMotion#(), epochrevolution&())
' Initialize beacon matrix and check age of satellite elements.
FOR i% = 1 TO numsats%
BeaconFrequency&(i%) = 100
IF epochyear%(i%) < fnyear% - 1 THEN
PRINT "Warning: Elements for satellite "; trim(name$(i%)); " are very old."
PRINT "You should update your "; tlefile$; " file."
END IF
NEXT
' Initialize the ground station data
CALL readgrnd
' ****** DRAW AND STORE SATELLITE SYMBOLS ******
CLS
SCREEN 2
LINE (4, 0)-(4, 4)
LINE (0, 2)-(8, 2)
GET (0, 0)-(8, 4), sat%
PUT (0, 0), sat%
CLS
LINE (4, 1)-(4, 3)
LINE (3, 2)-(5, 2)
GET (0, 0)-(8, 4), pkt%
PUT (0, 0), pkt%
' ***** Initializations done. *****
50
' ****** DERIVED CONSTANTS ******
' These values depend on the ground station data.
' C$=GROUND STATION CALL+LOCATION STRING
c$ = trim$(gsid$) + " " + trim$(gsloc$)
pi2# = 2 * pi#
SinGSLat# = SIN(fnradians(gslat!)) ' sin/cos of ground station latitude
CosGSLat# = COS(fnradians(gslat!))
SinGSLong# = SIN(fnradians(-gslong!)) ' sin/cos of ground station longitude
CosGSLong# = COS(fnradians(gslong!))
R9 = r0# * (1 - (f# / 2) + (f# / 2) * COS(2 * fnradians(gslat!))) + gsheight% / 1000
L8 = ATN((1 - f#) ^ 2 * SinGSLat# / CosGSLat#)
GSzcoord = R9 * SIN(L8)
GSxcoord = R9 * COS(L8) * CosGSLong#
GSycoord = R9 * COS(L8) * SinGSLong#
80
CLS
SCREEN 0
KEY(9) OFF
KEY(10) OFF
PRINT " USAT - QBASIC Version"
PRINT
PRINT "======================================================================"
PRINT
PRINT " SELECT ONE OF THE FOLLOWING OPTIONS:"
PRINT
PRINT " (P) ORBITAL PREDICTION PROGRAM"
PRINT
PRINT " (R) REALTIME TRACKING AND HIGH RESOLUTION SCREEN"
PRINT
PRINT " (S) Display ELEMENTS OF SATELLITES"
PRINT
PRINT " (G) CHANGE OR ENTER GROUNDSTATION DATA"
PRINT
PRINT " (D) RETURN TO DOS"
PRINT
PRINT
PRINT "ENTER SELECTION (P,R,C,G,D)--> "
CALL getkey(k$)
ON INSTR("PRSGD", k$) GOTO 4820, 2090, 240, 1620, 9000
BEEP
GOTO 50
240
' ****** SATFILE.BAS - VERSION 1.0, ISSUE 1.0 - HSIMODIF.1/25/85 ******
CLS
DO
PRINT "Elements of the following satellites are available:"
PRINT
FOR j% = 1 TO numsats%
PRINT name$(j%),
NEXT
PRINT
' ****** FIND SATELLITE RECORD ******
DO
PRINT : INPUT "Which satellite"; v1$: IF v1$ = "" THEN 80
v1$ = UCASE$(trim$(v1$))
FOR j% = 1 TO numsats%
IF UCASE$(trim$(name$(j%))) = v1$ THEN
CLS
' ****** DISPLAY SATELLITE RECORD ******
PRINT "Satellite = "; name$(j%)
PRINT "Epoch year = "; epochyear%(j%)
PRINT "Epoch day = "; epochjulianday#(j%)
PRINT "Inclination = "; inclination#(j%)
PRINT "R.A.A.N. = "; epochRAAN#(j%)
PRINT "Eccentricity = "; eccentricity#(j%)
PRINT "Arg. of perigee = "; epochArgOfPerigee#(j%)
PRINT "Mean anomaly = "; epochMA#(j%)
PRINT "Mean motion = "; MeanMotion#(j%)
PRINT "Epoch orbit no. = "; epochrevolution&(j%)
PRINT "Beacon freq. = "; BeaconFrequency&(j%)
PRINT
EXIT DO
END IF
NEXT
PRINT "Satellite not found."
LOOP
LOOP
1620
' ******* Groundsation data change v.1.0 OE1HSI jan.-1985**********
CLS
PRINT "CURRENT GROUND STATION DATA"
PRINT
GOSUB 2080
DO
PRINT "Do you want to CHANGE this DATA ? (Y/N)"
CALL getkey(k$)
ON INSTR("YN", k$) GOTO 1780, 2060
BEEP
LOOP
1710
PRINT
PRINT
GOSUB 2080
DO
PRINT "Do you want to make further changes ? (Y/N) "
CALL getkey(k$)
ON INSTR("YN", k$) GOTO 1780, 2050
BEEP
LOOP
1780
PRINT
PRINT "ENTER NEW DATA OR <RETURN> FOR UNCHANGED DATA":
PRINT
PRINT "Ground Station Identifier : ";
LINE INPUT u$
IF trim$(u$) <> "" THEN gsid$ = trim$(u$)
PRINT "Location of station : ";
LINE INPUT u$
IF trim$(u$) <> "" THEN gsloc$ = trim$(u$)
PRINT
PRINT "Longitude WEST of Greenwich (max +360) or East of Greenw. entered as -0 to -180"
PRINT
INPUT "Enter (with decimals) : ", u$
IF trim$(u$) <> "" THEN gslong! = VAL(u$)
IF gslong! < 0 THEN gslong! = 360 + gslong!
PRINT
PRINT "LATITUDE NORTH of Equator + (max 90) SOUTH of Equator - (max 90)"
PRINT
INPUT "ENTER (With decimals) : ", u$
IF trim$(u$) <> "" THEN gslat! = VAL(u$)
PRINT
INPUT "Station height above sea level in meters : ", u$
IF trim$(u$) <> "" THEN gsheight% = VAL(u$)
PRINT
PRINT "Correction from local computer time to UTC"
PRINT
INPUT "Enter (in hours with decimals) : ", u$
IF trim$(u$) <> "" THEN tzcor! = VAL(u$)
GOTO 1710
2050
PRINT
PRINT "DATA SAVED AS USATGRND.DAT"
OPEN "USATGRND.DAT" FOR OUTPUT AS #1
PRINT #1, gsid$
PRINT #1, gsloc$
PRINT #1, gslong!, gslat!, gsheight%, tzcor!
CLOSE #1
GOTO 50
2060
PRINT
PRINT "DATA NOT CHANGED"
GOTO 50
2080
PRINT "Ground Station Identifier : "; gsid$
PRINT "Location name : "; gsloc$
PRINT USING "West longitude (degrees) : +###.##"; gslong!
PRINT USING "Latitude (degrees) : +##.##"; gslat!
PRINT USING "Height above Mean Sea Level (meters) : #####"; gsheight%
PRINT USING "Local time correction to UTC (hours) : +##.##"; tzcor!
PRINT
RETURN
'**** END PROGRAM GROUNDSTATION DATA CHANGE/STORAGE OE1HSI JAN. 1985 ****
2090
' ****** ORBITS2 - VERSION 1.0, ISSUE 1.2 -11/1/83 ******
KEY OFF
SCREEN 2, 0
CLS
' Initialize D$ and T$ from date$ and time$
CALL timezone(DATE$, TIME$, d$, t$, tzcor!)
' ****** SET UP KEY TRAPPING ******
ON KEY(9) GOSUB 4760
KEY(9) STOP
ON KEY(10) GOSUB 4790
KEY(10) OFF
flg9 = 0
flg10 = 0
GOSUB 4290
' ****** ORBIT DETERMINATION LOOP STARTS HERE ******
q$ = ""
DO UNTIL q$ = CHR$(27)
FOR j% = 1 TO numsats%
q$ = UCASE$(INKEY$)
IF q$ = CHR$(27) THEN EXIT FOR
CALL timezone(DATE$, TIME$, d$, t$, tzcor!)
t = INT(30.55 * (fnmonth% + 2)) - 2 * (INT(.1 * (fnmonth% + 7))) - 91
IF fnmonth% > 2 AND fnleap%(fnyear%) THEN t = t + 1
IF epochyear%(j%) = fnyear% - 1 THEN
t = t + 365
IF fnleap%(k%) THEN t = t + 1
END IF
t = t + fnday% + fnhour% / 24# + fnminute% / 1440# + fnsecond% / 86400#
GOSUB 6780
IF flg9 THEN CALL plotloc(longitude!, latitude!, plotx%, ploty%): GOSUB 3890
GOSUB 4040
NEXT
LOOP
GOTO 50
' ****** PUT SATELLITE ON SCREEN ROUTINE ******
3890
GET (plotx% - 4, ploty% - 2)-(plotx% + 4, ploty% + 2), kep%
PUT (plotx% - 4, ploty% - 2), sat%, PRESET
PUT (plotx% - 4, ploty% - 2), sat%
PUT (plotx% - 4, ploty% - 2), kep%, PSET
PUT (plotx% - 4, ploty% - 2), sat%
IF FLG0 <> 0 THEN
PUT (satl(j%, 1), satl(j%, 2)), sat%
PUT (satl(j%, 1), satl(j%, 2)), pkt%, OR
END IF
satl(j%, 1) = plotx% - 4
satl(j%, 2) = ploty% - 2
IF j% = numsats% THEN FLG0 = 1
RETURN
' ****** PRINT SATELLITE DETAILS ROUTINE ******
4040
KEY(9) ON
KEY(10) ON
4050
KEY(10) STOP
KEY(9) STOP
IF FLGK = 1 THEN GOSUB 4270
IF flg9 = 0 THEN
LOCATE 2, 49
PRINT d$
LOCATE 2, 63
PRINT t$
IF elevation! >= 0 THEN
IF elevation! > 0 AND elevation! < 1 THEN BEEP
END IF
IF j% + 6 > 23 GOTO 4230
LOCATE j% + 6, 15
PRINT SPACE$(50)
LOCATE j% + 6, 8 - (LEN(trim$(name$(j%))) / 2)
PRINT trim$(name$(j%))
LOCATE j% + 6, 15
ELSE
IF flg10 = 1 THEN GOSUB 4540
IF flg10 <> 0 THEN
IF UCASE$(trim$(name$(j%))) <> u$ THEN RETURN
LOCATE 25, 69
PRINT "SELECTED";
END IF
LOCATE 25, 1
PRINT SPACE$(68);
LOCATE 25, 8 - (LEN(trim$(name$(j%))) / 2)
PRINT trim$(name$(j%));
LOCATE 25, 15
END IF
PRINT USING "### ### ##### ##### ###.# ###.# ######"; azimuth!; elevation!; range#; (r# - r0#); latitude!; longitude!; revolution&;
4230
IF flg9 <> 0 THEN LOCATE 20, 48: PRINT t$;
RETURN
' ****** SET UP SCREEN DISPLAY ROUTINE ******
4270
CLS
FLG0 = 0
FLGK = 0
4290
IF flg9 = 1 THEN
SCREEN 2, 0
DEF SEG = &HB800
BLOAD "USATCGA.MAP", 0
DEF SEG
CALL plotloc(gslong!, gslat!, plotx%, ploty%)
CIRCLE (plotx%, ploty%), 2
LOCATE , , 1, 12, 13
L1 = 22
L2 = 23
L3 = 24
LOCATE 20, 3
PRINT "Data for Groundstation At Time= UTC On: "; d$
LOCATE 20, 26
PRINT gsid$
ELSE
' FLG9=0
SCREEN 0, 1
LOCATE 1, 40 - LEN(c$) / 2, 0
PRINT c$
LOCATE 2, 5
PRINT "Real-time satellite tracking coordinates on at UTC"
LOCATE 25, 3
PRINT "F9 TOGGLES THE GRAPH/TABLE F10 TO SELECT SINGLE SAT IN GRAPH ESC TO END";
L1 = 4
L2 = 5
L3 = 6
END IF
LOCATE L1, 3
PRINT " NAME OR AZ EL RANGE HEIGHT LAT. LONG. ORBIT"
LOCATE L2, 3
PRINT "DESIGNATOR DEG DEG KM KM DEG DEG NO."
LOCATE L3, 3
PRINT "---------- --- --- ----- ------ ----- ----- ------";
RETURN
' ****** SELECT SINGLE SATELLITE ROUTINE ******
4540
LOCATE 25, 1
PRINT SPACE$(79);
LOCATE 25, 1
INPUT ; "WHICH SATELLITE? (CR FOR ALL)"; i$
i$ = UCASE$(trim$(i$))
IF i$ = "" THEN
flg10 = 0
ELSE
u$ = i$
flg10 = 2
END IF
LOCATE 25, 1
PRINT SPACE$(79);
RETURN
4760 ' toggle FLG9
flg9 = -(flg9 = 0)
FLGK = 1
RETURN 4050
4790 flg10 = flg9
RETURN 4050
4820
'****** ORBIT2 - VERSION 2.0, ISSUE 1.0/HSI - 17/01/85 *****
KEY OFF
SCREEN 0
CLS
' ****** HOUSEKEEPING ITEMS ******
C8$ = CHR$(10) + CHR$(10) + CHR$(10) + CHR$(10)
DO
PRINT "Enter start date (mm-dd-yy): ";
u$ = gets$
year% = VAL(RIGHT$(u$, 2))
month% = VAL(LEFT$(u$, 2))
day% = VAL(MID$(u$, 4, 2))
IF month% >= 1 AND month% <= 12 AND day% >= 1 AND day% <= 31 AND year% >= 91 AND MID$(u$, 3, 1) = "-" AND MID$(u$, 6, 1) = "-" THEN EXIT DO
LOOP
YY = year% + 1900
t$ = fnt$(year%) + "/" + fnt$(month%) + "/" + fnt$(day%) + " at "
D8 = day% + INT(30.55 * (month% + 2)) - 2 * (INT(.1 * (month% + 7))) - 91
IF month% > 2 THEN IF fnleap%(year%) THEN D8 = D8 + 1
PRINT " Day #"; D8
PRINT
PRINT "Enter start time (hh:mm): ";
DO
u$ = gets$
h% = VAL(LEFT$(u$, 2))
m% = VAL(RIGHT$(u$, 2))
IF h% >= 0 AND h% <= 23 AND m% >= 0 AND m% <= 59 AND MID$(u$, 3, 1) = ":" THEN EXIT DO
PRINT "Valid format is: hh:mm"
PRINT "Re-enter: ";
LOOP
T7 = D8 + h% / 24# + m% / 1440#
t$ = t$ + fnt$(h%) + fnt$(m%) + ":00 H"
PRINT "Enter duration in hours and minutes (hh:mm): ";
DO
u$ = gets$
h% = VAL(LEFT$(u$, 2))
m% = VAL(RIGHT$(u$, 2))
IF h% >= 0 AND m% >= 0 AND m% <= 59 AND MID$(u$, 3, 1) = ":" THEN EXIT DO
PRINT "Valid format is: hh:mm"
PRINT "Re-enter: ";
LOOP
T8 = T7 + h% / 24# + m% / 1440#
PRINT USING " From ###.####### to ###.#######"; T7; T8
DO
PRINT "Enter time increment (minutes): ";
INPUT m%
IF m% > 0 THEN EXIT DO
LOOP
T9 = m% / 1440#
PRINT
INPUT "MINIMUM ELEVATION ? (DEFAULT 0) Deg. = ", E8
DO
PRINT
INPUT "Output to Printer (P) or Screen (S) ?-->", P$
P$ = UCASE$(P$)
IF (P$ = "P" OR P$ = "S") THEN
EXIT DO
ELSE BEEP
END IF
LOOP
IF P$ = "P" THEN outfile$ = "lpt1:" ELSE outfile$ = "scrn:"
CLOSE #2
OPEN outfile$ FOR OUTPUT AS #2
IF outfile$ = "lpt1:" THEN
PRINT "Put the printer on-line and align the page,"
PRINT "then press any key when ready."
DO
getkey (u$)
IF u$ <> "" THEN EXIT DO
LOOP
END IF
' ****** SELECT SATELLITE FROM MENU ******
CLS
PRINT "SATELLITE SELECTION MENU"
PRINT
FOR j% = 1 TO numsats%
PRINT name$(j%),
NEXT
PRINT
DO
PRINT
INPUT "SELECT satellite >", y3$
y3$ = UCASE$(trim$(y3$))
FOR j% = 1 TO numsats%
IF UCASE$(trim$(name$(j%))) = y3$ THEN EXIT DO
NEXT
PRINT "Satellite not found."
LOOP
PRINT
PRINT "Doppler calculated for frequ. = "; BeaconFrequency&(j%); " MHz"
INPUT " Change frequency to: (0 for default) ", d
IF d <> 0 THEN BeaconFrequency&(j%) = d
PRINT #2,
PRINT #2, "Orbital ELEMENTS for "; name$(j%)
PRINT #2,
PRINT #2, "Reference epoch = "; epochyear%(j%); " +"; epochjulianday#(j%)
PRINT #2, "Starting epoch = "; year%; " +"; T7; " = "; t$
PRINT #2,
PRINT #2, "Parameter"; TAB(20); "Reference"; TAB(40); "Starting"
t = T7
IF epochyear%(j%) = year% - 1 THEN
t = t + 365
T8 = T8 + 365
IF fnleap%(k%) THEN t = t + 1: T8 = T8 + 1
END IF
' Calculate start time parameters
GOSUB 6780
PRINT #2, "Orbit Number "; TAB(20); epochrevolution&(j%); TAB(40); revolution&
PRINT #2, "Mean Anomaly "; TAB(20); epochMA#(j%); TAB(40); fndegrees(MeanAnomaly#)
PRINT #2, "Inclination "; TAB(20); inclination#(j%)
PRINT #2, "Eccentricity "; TAB(20); eccentricity#(j%)
PRINT #2, "Mean Motion "; TAB(20); MeanMotion#(j%)
PRINT #2, "S.M.A.,km "; TAB(20); SemiMajor#(j%)
PRINT #2, "Arg. Perigee "; TAB(20); epochArgOfPerigee#(j%); TAB(40); ArgOfPerigee#
PRINT #2, "R. A. A. N. "; TAB(20); epochRAAN#(j%); TAB(40); RAAN#
PRINT #2, "Freq.,MHz "; TAB(20); BeaconFrequency&(j%)
OldRevolution& = -1
OldJulianDay% = -1
'****** COMPUTATION LOOP ******
6400
t = t + T9
julianday% = INT(t)
IF julianday% <> OldJulianDay% THEN
OldJulianDay% = -1
OldRevolution& = -1
END IF
GOSUB 6780
IF elevation! >= E8 THEN
IF julianday% <> OldJulianDay% OR revolution& <> OldRevolution& THEN
IF julianday% <> OldJulianDay% THEN
GOSUB 6690
OldJulianDay% = julianday%
PRINT #2, " U.T.C. AZ EL DOPPLER RANGE HEIGHT LAT. LONG. PHASE"
PRINT #2, "HHMM:SS deg deg Hz km km deg deg <256>"
END IF
PRINT #2, TAB(21); "- - - ORBIT #"; revolution&; "- - -"
END IF
OldRevolution& = revolution&
T4 = t - julianday%
s4 = INT(T4 * 86400)
H4 = INT(s4 / 3600 + .000001)
M4 = INT((s4 - H4 * 3600) / 60 + .000001)
s4 = s4 - 3600 * H4 - 60 * M4
doppler# = -BeaconFrequency&(j%) * 1000000 * Vr / c#
PRINT #2, fnt$(H4) + fnt$(M4) + ":" + fnt$(s4);
PRINT #2, USING " ### ### #######"; azimuth!; elevation!; doppler#;
PRINT #2, USING " ##### ##### ###.# ###.# ###"; range#; (r# - r0#); latitude!; longitude!; phase%
END IF
IF t < T8 THEN GOTO 6400
PRINT #2,
PRINT "Do YOU have another INQUIRY ? (Y/N) "
PRINT
PRINT "Else you return to the MAIN MENU !"
PRINT
DO
CALL getkey(k$)
ON INSTR("YN", k$) GOTO 4820, 6670
BEEP
LOOP
6670
CLOSE #2
GOTO 80
'****** PAGE HEADER SUBROUTINE ******
6690
PRINT #2,
PRINT #2, c$; " Lat.="; gslat!; " W.Long.="; gslong!; " Ht.="; gsheight%;
PRINT #2, TAB(14); "- - - Minimum Elevation = "; E8; "Deg. - - -"
PRINT #2,
PRINT #2, TAB(14); "- - - DAY #"; julianday%; "- - - "; juliantodate$(julianday%, INT(YY)); " - - -"
PRINT #2,
RETURN
'****** ORBIT DETERMINATION AND UTILITY ROUTINES ******
' Generate position data within the satellite orbital plane
6780
' Calculate the Mean Anomaly (in radians) within the orbital plane
OrbitPos# = epochrevolution&(j%) + epochMA#(j%) / 360 + MeanMotion#(j%) * (t - epochjulianday#(j%))
revolution& = INT(OrbitPos#)
phase% = INT((OrbitPos# - revolution&) * 256)
MeanAnomaly# = (OrbitPos# - revolution&) * pi2#
' Newton-Raphson calculation of Eccentric Anomaly
EA# = MeanAnomaly# + eccentricity#(j%) * SIN(MeanAnomaly#) + .5 * (eccentricity#(j%) * eccentricity#(j%)) * SIN(2 * MeanAnomaly#)
DO
SinEA# = SIN(EA#)
CosEA# = COS(EA#)
slope# = 1 - eccentricity#(j%) * CosEA#
Yerror# = EA# - eccentricity#(j%) * SinEA# - MeanAnomaly#
IF ABS(Yerror#) < .000001 THEN EXIT DO ELSE EA# = EA# - Yerror# / slope#
LOOP
' Calculate the satellite position in the orbit plane
SemiMajor#(j%) = (g0# / (MeanMotion#(j%) * MeanMotion#(j%))) ^ (1 / 3)
ba2# = 1 - eccentricity#(j%) * eccentricity#(j%)
orbitX# = SemiMajor#(j%) * (CosEA# - eccentricity#(j%))
orbitY# = SemiMajor#(j%) * SQR(ba2#) * SinEA#
r# = SemiMajor#(j%) * slope#
' Calculate relative orientation of the satellite orbit
PrecessFactor# = 9.95 * ((r0# / SemiMajor#(j%)) ^ 3.5) / (ba2# * ba2#)
SinInclination# = SIN(fnradians(inclination#(j%)))
CosInclination# = COS(fnradians(inclination#(j%)))
RAAN# = epochRAAN#(j%) - (t - epochjulianday#(j%)) * PrecessFactor# * CosInclination#
ArgOfPerigee# = epochArgOfPerigee#(j%) + (t - epochjulianday#(j%)) * PrecessFactor# * (2.5 * (CosInclination# * CosInclination#) - .5)
' Generate the satellite orbital-plane to
' celestial coordinates translation matrix
SinRAAN# = SIN(fnradians(RAAN#))
CosRAAN# = COS(fnradians(RAAN#))
SinAoP# = SIN(fnradians(ArgOfPerigee#))
CosAoP# = COS(fnradians(ArgOfPerigee#))
cc#(1, 1) = (CosAoP# * CosRAAN#) - (SinAoP# * SinRAAN# * CosInclination#)
cc#(1, 2) = -(SinAoP# * CosRAAN#) - (CosAoP# * SinRAAN# * CosInclination#)
cc#(2, 1) = (CosAoP# * SinRAAN#) + (SinAoP# * CosRAAN# * CosInclination#)
cc#(2, 2) = -(SinAoP# * SinRAAN#) + (CosAoP# * CosRAAN# * CosInclination#)
cc#(3, 1) = (SinAoP# * SinInclination#)
cc#(3, 2) = (CosAoP# * SinInclination#)
' Calculate position in the celestial coordinate system (origin at geocenter)
celestialX# = orbitX# * cc#(1, 1) + orbitY# * cc#(1, 2)
celestialY# = orbitX# * cc#(2, 1) + orbitY# * cc#(2, 2)
celestialZ# = orbitX# * cc#(3, 1) + orbitY# * cc#(3, 2)
' Now take Earth's rotation and solar orbit into account
GreenwichDeclination# = t * g1# + fngmst#(1900 + epochyear%(j%))
GreenwichDeclination# = (GreenwichDeclination# - INT(GreenwichDeclination#)) * pi2#
SinGD# = SIN(GreenwichDeclination#)
CosGD# = COS(GreenwichDeclination#)
' Calculate sub-satellite position
SSPX# = (celestialX# * CosGD#) + (celestialY# * SinGD#)
SSPY# = -(celestialX# * SinGD#) + (celestialY# * CosGD#)
SSPZ# = celestialZ#
longitude! = 360 - fndegrees(arctan(SSPY#, SSPX#))
latitude! = fndegrees(fnarcsin(SSPZ# / r#))
XRange# = (SSPX# - GSxcoord)
YRange# = (SSPY# - GSycoord)
ZRange# = (SSPZ# - GSzcoord)
range# = SQR(XRange# * XRange# + YRange# * YRange# + ZRange# * ZRange#)
' Translate to ground station coordinate system
GSCSZ# = (XRange# * CosGSLong# * CosGSLat#) + (YRange# * SinGSLong# * CosGSLat#) + (ZRange# * SinGSLat#)
GSCSX# = -(XRange# * CosGSLong# * SinGSLat#) - (YRange# * SinGSLong# * SinGSLat#) + (ZRange# * CosGSLat#)
GSCSY# = (YRange# * CosGSLong#) - (XRange# * SinGSLong#)
elevation! = fndegrees(fnarcsin(GSCSZ# / range#))
azimuth! = fndegrees(arctan(GSCSY#, GSCSX#))
' Speed to/from ground station
IF OldT# <> t THEN Vr = ((range# - OldRange#) / (t - OldT#)) / 86400 ELSE Vr = -900000000
OldRange# = range#
OldT# = t
RETURN
9000 END
FUNCTION arctan (y, x)
' Modified arctangent function. Returns a value between 0 and 2 PI.
SHARED pi#, pi2#
IF x > 0 AND y >= 0 THEN
arctan = ATN(y / x)
ELSEIF x > 0 AND y < 0 THEN
arctan = ATN(y / x) + pi2#
ELSEIF x < 0 THEN
arctan = ATN(y / x) + pi#
ELSEIF y >= 0 THEN
arctan = pi# / 2
ELSE arctan = 1.5# * pi#
END IF
END FUNCTION
DEFSNG A-Z
' Converts a date string to the julian day of the year
FUNCTION datetojulian% (d$)
day% = VAL(MID$(d$, 4, 2))
SELECT CASE LEFT$(d$, 2)
CASE "01"
CASE "02"
day% = day% + 31
CASE "03"
day% = day% + 59
CASE "04"
day% = day% + 90
CASE "05"
day% = day% + 120
CASE "06"
day% = day% + 151
CASE "07"
day% = day% + 181
CASE "08"
day% = day% + 212
CASE "09"
day% = day% + 243
CASE "10"
day% = day% + 273
CASE "11"
day% = day% + 304
CASE "12"
day% = day% + 334
END SELECT
IF fnleap%(VAL(RIGHT$(d$, 4))) AND LEFT$(d$, 2) > "02" THEN
datetojulian% = day% + 1
ELSE
datetojulian% = day%
END IF
END FUNCTION
SUB getkey (k$)
k$ = ""
WHILE k$ = ""
k$ = UCASE$(INKEY$)
WEND
END SUB
' Gets an input line, trims leading and trailing blanks
' then converts to upper-case before returning.
FUNCTION gets$
LINE INPUT i$
gets$ = UCASE$(trim$(i$))
END FUNCTION
' Converts a julian day and year to a date string.
FUNCTION juliantodate$ (julianday%, year%)
day% = INT(julianday%)
IF day% = 60 THEN
IF fnleap%(year%) THEN tempdate$ = "02-29" ELSE tempdate$ = "03-01"
GOTO appendyear
ELSE IF fnleap%(year%) AND day% > 60 THEN day% = day% - 1
END IF
SELECT CASE day%
CASE 1 TO 31
tempdate$ = "01-" + MID$(STR$(day% + 100), 3)
CASE 32 TO 59
tempdate$ = "02-" + MID$(STR$(day% - 31 + 100), 3)
CASE 60 TO 90
tempdate$ = "03-" + MID$(STR$(day% - 59 + 100), 3)
CASE 91 TO 120
tempdate$ = "04-" + MID$(STR$(day% - 90 + 100), 3)
CASE 121 TO 151
tempdate$ = "05-" + MID$(STR$(day% - 120 + 100), 3)
CASE 152 TO 181
tempdate$ = "06-" + MID$(STR$(day% - 151 + 100), 3)
CASE 182 TO 212
tempdate$ = "07-" + MID$(STR$(day% - 181 + 100), 3)
CASE 213 TO 243
tempdate$ = "08-" + MID$(STR$(day% - 212 + 100), 3)
CASE 244 TO 273
tempdate$ = "09-" + MID$(STR$(day% - 243 + 100), 3)
CASE 274 TO 304
tempdate$ = "10-" + MID$(STR$(day% - 273 + 100), 3)
CASE 305 TO 334
tempdate$ = "11-" + MID$(STR$(day% - 304 + 100), 3)
CASE 335 TO 365
tempdate$ = "12-" + MID$(STR$(day% - 334 + 100), 3)
END SELECT
appendyear: juliantodate$ = tempdate$ + "-" + MID$(STR$(year%), 2)
END FUNCTION
SUB plotloc (longitude!, latitude!, xcoord%, ycoord%)
' Returns the screen x and y coordinates for a given map location
ycoord% = CINT(.7111 * (90 - latitude!) + 3)
IF longitude! >= 0 AND longitude! <= 270 THEN
xcoord% = CINT(477 - longitude! * 1.7444)
ELSE xcoord% = CINT(1105 - longitude! * 1.7444)
END IF
END SUB
SUB readgrnd
SHARED gsid$, gsloc$, gslong!, gslat!, gsheight%, tzcor!
' Format of USATGRND.DAT:
' gsid$ = Station Identifier
' gsloc$ = Station location (name)
' gslong! = Longitude in degrees
' gslat! = Latitude in degrees
' gsheight% = Height above sea level in meters
' tzcor! = Time zone correction to UTC
OPEN "USATGRND.DAT" FOR INPUT AS #1
LINE INPUT #1, gsid$
LINE INPUT #1, gsloc$
INPUT #1, gslong!, gslat!, gsheight%, tzcor!
CLOSE #1
END SUB
FUNCTION readtle% (tlefile$, name$(), epochyear%(), epochjulianday#(), inclination#(), epochRAAN#(), eccentricity#(), epochArgOfPerigee#(), epochMA#(), MeanMotion#(), epochrevolution&())
' ****** ESTABLISH SATELLITE ELEMENT MATRIX ******
' tlefile$ conforms to the NASA/NORAD three-line element set format.
OPEN tlefile$ FOR INPUT AS #1
FOR i% = 1 TO UBOUND(name$)
IF EOF(1) THEN EXIT FOR
LINE INPUT #1, name$(i%)
IF trim$(name$(i%)) = "" THEN EXIT FOR
LINE INPUT #1, y3$
LINE INPUT #1, t0$
epochyear%(i%) = VAL(MID$(y3$, 19, 2))
epochjulianday#(i%) = VAL(MID$(y3$, 21, 12))
inclination#(i%) = VAL(MID$(t0$, 9, 8))
epochRAAN#(i%) = VAL(MID$(t0$, 18, 8))
eccentricity#(i%) = VAL("." + MID$(t0$, 27, 7))
epochArgOfPerigee#(i%) = VAL(MID$(t0$, 35, 8))
epochMA#(i%) = VAL(MID$(t0$, 44, 8))
MeanMotion#(i%) = VAL(MID$(t0$, 53, 11))
epochrevolution&(i%) = VAL(MID$(t0$, 64, 5))
NEXT
CLOSE #1
readtle% = i% - 1
END FUNCTION
' Converts an input date and time string from one time zone to another.
SUB timezone (d0$, t0$, d1$, t1$, tzcor!)
julianday% = datetojulian(d0$)
year% = VAL(RIGHT$(d0$, 4))
hour! = VAL(LEFT$(t0$, 2)) + INT(tzcor!) + (VAL(MID$(t0$, 4, 2)) + tzcor! - INT(tzcor!)) / 60
SELECT CASE hour!
CASE IS < 0
hour! = hour! + 24
julianday% = julianday% - 1
CASE IS >= 24
hour! = hour! - 24
julianday% = julianday% + 1
END SELECT
t1$ = fnt$(INT(hour!)) + ":" + fnt$(INT(60 * (hour! - INT(hour!)))) + MID$(t0$, 6)
SELECT CASE julianday%
CASE 0
year% = year% - 1
IF fnleap%(year%) THEN julianday% = 366 ELSE julianday% = 365
CASE IS > 365
IF NOT fnleap%(year%) THEN julianday% = 1: year% = year% + 1
END SELECT
d1$ = juliantodate$(julianday%, year%)
END SUB
' Trims leading and trailing spaces from a string.
FUNCTION trim$ (i$)
trim$ = RTRIM$(LTRIM$(i$))
END FUNCTION